// Represents a square matrix
Matrix
	var
		rows
		columns
		list/values

	New(r, c = r, fill = 0)
		. = ..()
		ASSERT(r && c)
		src.rows = r
		src.columns = c
		src.values = new(src.rows * src.columns)
		for(var/i = 1 to src.values.len)
			src.SetValue(getRow(i), getCol(i), fill)
	proc
		getRow(index)
			return Math.Ceiling(index / src.columns)
		getCol(index)
			return index - (getRow(index) - 1) * columns
		getIndex(row, col)
			return col + (row - 1) * src.columns

		Add(Matrix/M)
			ASSERT(SameSize(M))
			var/Matrix/result = src.Copy()
			for(var/i = 1 to src.values.len)
				var
					row = getRow(i)
					col = getCol(i)
				result.SetValue(row, col, src.GetValue(row, col) + \
									M.GetValue(row, col))
			return result

		Copy()
			var/Matrix/M = new src.type(src.rows, src.columns)
			for(var/i = 1 to values.len)
				var
					row = getRow(i)
					col = getCol(i)
				M.SetValue(row, col,
					src.GetValue(row, col))
			return M


		Equals(Matrix/M)
			if(!SameSize(M))
				return FALSE
			for(var/i = 1 to values.len)
				var
					row = getRow(i)
					col = getCol(i)
				if(src.GetValue(row, col) !=  M.GetValue(row, col))
					return FALSE
			return TRUE

		IsSquare()
			return src.rows == src.columns

		Multiply(Matrix/M)
			var/Matrix/result
			if(isnum(M))
				result = new(rows, columns)
				for(var/i = 1 to values.len)
					var
						row = getRow(i)
						col = getCol(i)
					result.SetValue(row, col,
						src.GetValue(row, col) * M)
			else
				ASSERT(istype(M))
				ASSERT(src.columns == M.rows)
				result = new(src.rows, M.columns)
				for(var/row = 1 to result.rows)
					for(var/col = 1 to result.columns)
						var/sum = 0
						for(var/i = 1 to M.rows)
							sum += src.GetValue(row, i) * M.GetValue(i, col)
						result.SetValue(row, col, sum)
			return result

		Populate(...)
			ASSERT(args.len <= src.rows * src.columns)
			for(var/i = 1 to values.len)
				var
					row = getRow(i)
					col = getCol(i)
				src.SetValue(row, col, args[i])
				if(++i > args.len)
					return

		SameSize(Matrix/M)
			return M.rows == src.rows && M.columns == src.columns

		ToString()
			. = "\["
			for(var/row = 1 to rows)
				. += "\["
				for(var/col = 1 to columns)
					. += "[src.GetValue(row, col)]"
					if(col < columns)
						. += " "
				. += "]"
				if(row < rows)
					. += "\n"
			. += "]"

		Transpose()
			var/Matrix/result = new(columns, rows)
			for(var/i = 1 to values.len)
				var
					row = getRow(i)
					col = getCol(i)
				result.SetValue(col, row, src.GetValue(row, col))
			return result

		GetValue(row, col)
			ASSERT(row <= rows && col <= columns)
			return src.values[getIndex(row, col)]

		SetValue(row, col, value)
			ASSERT(row <= rows && col <= columns && isnum(value))
			src.values[getIndex(row, col)] = value

	// Square matrices only.
	// It's easier to have them under the same type as normal matrices
	// so that if an operation happens to return a square matrix,
	// it does not need to be copied over to a new object in order to perform these
	// operations.
		Adjugate()
			ASSERT(src.IsSquare())
			var/Matrix/M = src.CofactorMatrix()
			return M.Transpose()

		Cofactor(row, col)
			ASSERT(src.IsSquare())
			var/m = src.Minor(row, col)
			if(Math.IsOdd(row + col)) m *= -1
			return m

		CofactorMatrix()
			ASSERT(src.IsSquare())
			var/Matrix/M = new(src.rows)
			for(var/i = 1 to values.len)
				var
					row = getRow(i)
					col = getCol(i)
				M.SetValue(row, col, src.Cofactor(row, col))
			return M

		Determinant()
			ASSERT(src.IsSquare())
			if(src.rows == 1)
				return src.GetValue(1, 1)
			var/sum = 0
			for(var/col = 1 to columns)
				sum += src.GetValue(1, col) * src.Cofactor(1, col)
			return sum

		Inverse()
			ASSERT(src.IsSquare())
			var/Matrix/M = src.Adjugate()
			var/d = src.Determinant()
			if(!d) return null
			return M.Multiply(1 / d)

		Minor(row, col)
			ASSERT(src.IsSquare())
			var/Matrix/M = new(src.rows - 1)
			for(var/i = 1 to rows)
				if(i == row)
					continue
				for(var/j = 1 to columns)
					if(j == col)
						continue
					var
						mi 	= i < row ? i : i - 1
						mj 	= j < col ? j : j - 1
					M.SetValue(mi, mj, src.GetValue(i, j))
			return M.Determinant()

// Matrix helper methods
var/_matrix/Matrix = new()
_matrix
	proc
		Identity(size)
			var/Matrix/M = new(size)
			for(var/i = 1 to size)
				M.SetValue(i, i, 1)
			return M


//Optimized implementation of a 2x2 matrix.
Matrix22
	parent_type = /Matrix
	var
		M11; M12
		M21; M22

	New()
		. = ..(2)

	Determinant()
		return M11 * M22 - M21 * M12

	Inverse()
		var
			Matrix22/M  = src.Multiply(1 / src.Determinant())
			tM11		= M.M11
		M.M11 = M.M22
		M.M22 = tM11
		M.M12 = -M.M12
		M.M21 = -M.M21
		return M

	GetValue(row, col)
		ASSERT(row <= rows && col <= columns)
		return src.vars["M[row][col]"]

	SetValue(row, col, value)
		ASSERT(row <= rows && col <= columns && isnum(value))
		src.vars["M[row][col]"] = value

/*
Matrix33
	parent_type = /Matrix
	size		= 3

	var
		M11; M12; M13
		M21; M22; M23
		M31; M32; M33

	Determinant()
		return M11 * M22 * M33 + M12 * M23 * M31 + M13 * M21 * M32 \
			 - M31 * M22 * M13 - M32 * M23 * M11 - M33 * M21 * M12
*/